Fast paired method of a 1-D cyclic convolution realization

ABSTRACT

Systems and methods are described for a fast paired method of 1-D cyclic convolution. A method includes calculating a paired transform of a signal, grouping components of the paired transform to form a plurality of splitting-signals, shifting the plurality of splitting signals, multiplying the plurality of splitting signals by a plurality of corresponding Fourier transforms, and calculating an inverse paired transform of the plurality of splitting signals.

BACKGROUND OF THE INVENTION

1. Field of the Invention

The invention relates generally to the field of digital signalprocessing. More particularly, the invention relates to digitalfiltering for digital signal processing.

2. Discussion of the Related Art

Digital signal processing is an important aspect of digitalcommunications. It may be used in communication applications such assatellite communication systems, satellite ranging systems (GPS, WAAS,GLONAS), radar systems, radio-locators for the military, and speechprocessing (such as processing a speech waveforms, when the input signalis of infinite duration), and restoration of a signal from noisy input.

Three important tools in signal processing are FIR filters, optimalWiener filters, and correlation functions. Linear convolution is anothertool that is realized in many processors, including industriallyavailable processors such as: Universal processors and mobile devices'processors from Texas Instruments (Dallas, Tex.) and ST Microelectronics(Dallas, Tex.).

FIR Filters

Finite-length impulse response (FIR) filters are widely used. FIRfilters may be important elements of most digital systems for thepreviously mentioned applications.

An 1-D FIR filter may be described by the linear convolution typeequation

$\begin{matrix}{{{y(n)} = {\sum\limits_{k = 0}^{M - 1}\;{{x( {k - n} )}{h(k)}}}},\mspace{31mu}{n = {0:( {N - 1} )}},} & (1)\end{matrix}$where x(n) is an input signal of length L and h(k) is the impulseresponse characteristics of the filter that has a length M≧L. IfN≧L+M−1, the cyclic convolution and linear convolution are identical.

As is known, the efficient realization of FIR filters is not provided bycalculation in the time domain, especially for long data signals. Thefrequency-domain approach based on the discrete Fourier transform is acomputational procedure that serves as an alternative to time-domainconvolution, and it is computationally more efficient than time-domainconvolution due to the existence of efficient algorithms for computingthe DFT. These algorithms are called fast Fourier transform (FFT)algorithms.

Optimal Wiener Filters

A Wiener filter is defined as a filter that provides optimal estimate oforiginal discrete signal x(n) from an observed noisy signal

$\begin{matrix}{{{y(n)} = {{\sum\limits_{k = 0}^{N - 1}\;{{x( {k - n} )}{h(k)}}} + {n(n)}}},\mspace{31mu}{n = {0:( {N - 1} )}},} & (2)\end{matrix}$where n(n) is noise. The optimality is with respect to mean-square errorand the estimate is defined as the linear cyclic convolution

$\begin{matrix}{{{\hat{x}(n)} = {\sum\limits_{k = 0}^{N - 1}\;{{w(k)}{y( {k - n} )}}}},\mspace{31mu}{n = {0:( {N - 1} )}},} & (3)\end{matrix}$where w(k) is the impulse response function of the Wiener filter. Thefiltering in equation (3) is implemented in the Fourier domain,{circumflex over (X)}(ω)=W(ω)Y(ω), by the filter with transfer function

$\begin{matrix}{{{W(\omega)} = \frac{\overset{\_}{H}(\omega)}{{{H(\omega)}}^{2} + {\phi_{n/x}(\omega)}}},\mspace{31mu}{\omega = {0:( {N - 1} )}},} & (4)\end{matrix}$where H(ω) is the DFT of h(k) and φ_(n/x) (ω), the noise-signal ratio.{circumflex over (X)}(ω) and Y(ω) are the Fourier transforms of{circumflex over (x)}(n) and y(n) respectively.

The realization of optimal filtering requires three DFTs and Nmultiplications to calculate {circumflex over (X)}(ω). Therefore, thetotal number of multiplications required to perform optimal filtrationby the DFT method can be estimated as3[N/2(log₂ N−3)+2]+N+3N=3N/2(log₂ N−3)+4N+6.where 3N operations are assumed for calculation of the optimal filter inequation (4).Cross-Correlation

Cross-correlation between two discrete sequences or signals x(n) andh(n) is defined as

$\begin{matrix}{{r(n)} = {\frac{1}{N}{\sum\limits_{k = 0}^{N - 1}\;{{x(k)}{h( {k + n} )}}}}} & (5)\end{matrix}$

The data sequences may be either correlated or convolved simply byreversing the order of one of the data sequences. Thus, the correlationand convolution can be calculated by the same program simply byreversing one of the sequences.

The correlation computation may be speeded by using the DFT approachbecause in the Fourier domain the correlation takes the following formR(ω)=X(ω)Ĥ(ω), ω=0:(N−1).  (6)

The correlation realization by the FFT requires the same number ofoperations as for the linear convolution. For longer data sequences, thefast DFT method may be used to achieve fast correlation.

FFT Method of Convolution Realization

The DFT provides a discrete frequency representation of thefinite-duration sequence (signal) in frequency domain and it reduces thecyclic convolution of two sequences to the simple operation of themultiplication of their Fourier transformsY(p)=X(p)H(p)  (7)where Y(p), X(p), and H(p) are the discrete Fourier transforms of y(n),x(n), and h(k), respectively.

The discrete Fourier transform, X(p), of the signal x(t), of length N,is multiplied by the frequency characteristics, H(p), of a linearfilter. And then the inverse discrete Fourier transform, F⁻¹, is used,which results in the filtered signal, Y(t). A diagram of the Fouriermethod is shown in blocks 100-105 depicted in FIG. 1A.r(t)→F[r](p)=X=(p)→X(p)H(p)=Y(p)→F ⁻¹ [Y(p)](t)=y(t)

The DFT may be performed twice (the DFT performance is not consideredfor the H(p)). This approach uses N values of the frequencycharacteristics and requires N operations of complex multiplication inthe frequency domain and N−1 additions. Additionally, direct and inverseDFTs require N/2 (log₂ N−3)+2 complex multiplications. Therefore, thetraditional FFT method of the linear filtration requiresM _(N)=2[N/2(log₂ N−3)+2]+N=N(log₂ N−2)+4.  (8)complex multiplications, and aboutA _(N)=2N(log₂ N−2)  (9)operations of additions. Here, the case of most interest in practice isconsidered, when the length of the signal is N=2^(r), r>1.

If one assumes that the complex multiplication is performed by threereal multiplications and three additions, then the numbers of realoperations of multiplications and additions can be estimated as3M _(N)=3N(log₂ N−2)+9, 3A _(N)=6N(log₂ N−2).  (10)The FFT-based method of realization of the convolution is superior froma computational point of view.Other Methods

Some disadvantages of applying the DFT are in the use of trigonometricfunctions cosine and sine functions, and other complex arithmeticfunctions. To avoid such disadvantages, different algebraic methods havebeen proposed. Among them are the overlap-saved and overlap-add methodsand the nesting of polynomials. These algorithms represent 1-Dconvolution as multi-dimensional convolutions which are reduced tocomputation of small length convolutions and polynomial multiplications.These algorithms require fewer arithmetic operations than the FFTapproach for small sequence length N<200. This is due to the fact thatthe number of multiplications per sample for these algorithms equals(3/2)^(r) for N=2^(r), and this number grows exponentially faster thanthe linear estimate (3r−2) for the DFT. Therefore, these algorithms donot provide the fast performance of the FIR filter for long data.

The requirements for a fast and simple convolution process for longsignals have not been fully met. What is needed is a solution thataddresses both of these requirements.

SUMMARY OF THE INVENTION

There is a need for the following embodiments. Of course, the inventionis not limited to these embodiments.

According to an aspect of the invention, a method for cyclic convolutioncomprises calculating a paired transform of a signal, groupingcomponents of the paired transform to form a plurality ofsplitting-signals, shifting the plurality of splitting signals,multiplying the plurality of splitting signals by a plurality ofcorresponding Fourier transforms, and calculating an inverse pairedtransform of the plurality of splitting signals.

BRIEF DESCRIPTION OF THE DRAWINGS

The drawings accompanying and forming part of this specification areincluded to depict certain aspects of the invention. The invention maybe better understood by reference to one or more of these drawings incombination with the description presented herein. It should be notedthat the features illustrated in the drawings are not necessarily drawnto scale.

FIG. 1A is a flow diagram showing the traditional method of signalfiltration by a Discrete Fourier Transform (DFT).

FIG. 1B is a flow diagram showing a paired-method of signal filtrationby the discrete paired transform (DPT), in accordance with an embodimentof the present disclosure.

FIG. 2 illustrates the basis functions of the 2, 4, and 8-point pairedtransforms, an aspect of the present disclosure.

FIG. 3 illustrates the splitting of a 256-point signal by a pairedtransform, in accordance with an embodiment of the present disclosure.

FIG. 4 illustrates angles |φ_(p)| and φ_(p) on a circle.

FIG. 5 shows a 256-point linear signal filtration by a paired transform,in accordance with an embodiment of the present disclosure.

FIG. 6 shows Fourier characteristics of a linear filter.

DETAILED DESCRIPTION

Linear convolution is the core of digital signal and image processing. Anew fast method for implementing one-dimensional cyclic (circular)convolution is presented. The method is based on the pairedrepresentation of a 1-D signal with respect to the Fourier transform. Asignal is transformed uniquely into a set of short signals, by means ofwhich the 1-D DFT of the signal is split into a set of short 1-D DFTs.

The paired transformation allows one to perform a linear filtration of asignal by means of an incomplete number of the frequency characteristicsof the filter at specific points. As a result, the complexity of theimplementation decreases. For example, the linear filtration of thesignal of length N=2^(r), r>1, can be performed by using only log₂ N+1values of frequency characteristics of the filters and N realmultiplications, instead of N log₂ N complex multiplications in thetraditional method.

Linear convolution of two sequences x(n) and h(n) of length N is definedby equation

$\begin{matrix}{{{y(n)} = {\sum\limits_{k = 0}^{N - 1}{{x( {k - n} )}{h(k)}}}},\mspace{31mu}{n = {0:( {N - 1} )}},} & (11)\end{matrix}$where indices are taken modulo N. The direct calculation of theconvolution requires N² multiplications.Fast Method of Paired Transforms

This disclosure demonstrates an approach that is based on the pairedrepresentation of signals, which is a signal transformation to a set ofshort signals. This transformation, which is called the pairedtransformation, allows for the performance of linear convolution(filtration) of the signal by performing convolutions (filtration) ofshort signals with fewer operations than is required by the method basedon DFT.

The proposed fast paired method of the convolution may be applied atleast for:

1. Fast realization of noncyclic linear convolutions, or FIR filters;

2. Fast realization of the optimal Wiener filtration of signals; and

3. Fast realization of the cyclic correlation.

Instead of filtering signals in the frequency domain, linear convolution(filtration) is implemented by processing the short signals. Such methodrequires only N operations of real multiplication and uses only a few ofthe frequency characteristics of a linear filter at specific points. Asa result, the paired method of linear convolution realization uses log₂N operations less than the FFT method. For instance, the linearconvolution of the signal of length 1024 by the paired transforms uses10 times fewer operations than the FFT method requires.

A diagram of the realization of the linear convolution by the pairedtransforms, an embodiment of the present invention, is given in FIG. 1B.As shown in blocks 100-111 of FIG. 1B, multiplication steps are takenonly at a few selected points of the filter Y(p), instead of multiplyingthe signal by all values of the filter, as shown in FIG. 1A. To describethe proposed method, the definition and main properties of the unitarypaired transforms that are used in the method are given.

Unitary Paired Transform

The one-dimensional unitary discrete paired transforms (DFT's) isdescribed in the following way. Consider the case of most interest, whenthe length of signals is N=2^(r), r>1. Let p,tεX_(N)={0,1, . . . , N−1},and let χ_(p,t) be the binary function

$\begin{matrix}{{x_{p,t}(n)} = \{ {{\begin{matrix}{1,} & {{{{if}\mspace{14mu}{np}} = {t\mspace{14mu}{mod}\mspace{14mu} N}},} \\{0,} & {{otherwise},}\end{matrix}\mspace{14mu} n} = {0:{( {N - 1} ).}}} } & (12)\end{matrix}$

Definition: Given a sample pεX_(N) and integer tε[0, N/2], the functionχ′_(p,t)(n)=X_(p,t)(n)−X_(p,t+N/2)(n)  (13)is called the 2-paired, or shortly, the paired function.

The totality of the paired function{χ′_(2n,2nt) ;n=0:(r−1),t=0:(2^(r−n−1)−1),1}  (14)is the complete and orthogonal set of functions. The unitary transformχ′ composed by this complete set of functions is called the pairedtransform. For example, the matrix of the 8-point paired transform χ′₈takes the form

$X_{8}^{1} = {\begin{bmatrix}\lbrack x_{1,0}^{\prime} \rbrack \\\lbrack x_{1,1}^{\prime} \rbrack \\\lbrack x_{1,2}^{\prime} \rbrack \\\lbrack x_{1,3}^{\prime} \rbrack \\\lbrack x_{2,4}^{\prime} \rbrack \\\lbrack x_{2,2}^{\prime} \rbrack \\\lbrack x_{4,0}^{\prime} \rbrack \\\lbrack x_{n,0}^{\prime} \rbrack\end{bmatrix} = \begin{bmatrix}1 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 \\0 & 1 & 0 & 0 & 0 & {- 1} & 0 & 0 \\0 & 0 & 1 & 0 & 0 & 0 & {- 1} & 0 \\0 & 0 & 0 & 1 & 0 & 0 & 0 & {- 1} \\1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} & 0 \\0 & 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} \\1 & {- 1} & 1 & {- 1} & 1 & {- 1} & 1 & {- 1} \\1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\end{bmatrix}}$Splitting of the DFT

Here, the one-dimensional N-point Fourier transform of a discrete signalf={f_(n)}_(n=0:(N−1)) is considered:

$\begin{matrix}{{F_{p} = {( {\mathcal{F}_{N} \circ f} )_{p} = {\sum\limits_{n = 0}^{N - 1}\;{f_{n}W^{np}}}}},{p = {0:( {N - 1} )}},} & (15)\end{matrix}$where W=W_(N)=exp(−j2π/N).

By means of N/2^(k)-point paired transforms χ′_(N/2) _(k) , k=0:(r−1),the matrix of the Fourier transform is decomposed into the products ofthe weakly filled binary matrices. The paired transform is fast andrequires 3N−2 additions. The decomposition of the N-point DFT by thepaired transforms requires 2^(r−1)(r−3)+2 multiplications.

For each point pεX_(N), the following equality holds

$\begin{matrix}{{F_{\overset{\_}{{({{2m} + 1})}p}} = {\sum\limits_{t = 0}^{{N/2} - 1}\;{( {f_{p,t}^{\prime}W^{t}} )W_{N/2}^{mt}}}},\mspace{34mu}{m = {0:{( {{N/2} - 1} ).}}}} & (16)\end{matrix}$

The sequence {f′_(p,0), f′_(p,1), . . . , f′_(p, N/2−1)} determinescompletely the N-point DFT of the signal f at points of the set T′_(p)of X_(N),T′ _(p)={ (2m+1)p ; m=0:(N/2−1)}where l=l mod N, for an integer l:

Such a sequence may be considered a splitting-signal of the signal f,that corresponds to the set T′_(p). Consider the following partition σ′of X_(N) composed of the sets T′_(p)σ′=σ′_(N)=((T′ ₂ _(n) ;n=0:r−1),{0})  (17)

The splitting-signal that corresponds to the first set T′₁ of thepartition or σ′f′_(T′) ₁ ={f′_(1,0), f′_(1,1), . . . , f′_(1,N/2−1)}has length N/2. For the next splitting-signal f_(T′) ₂ , all componentsf′_(2,t)=0 with odd t, and equation (16) represents the N/4-point DFT(note, that the set T′₂ contains also N/2 points). Therefore,noninformative components f′_(2,2k+1) can be discarded, and f′_(T) ₂ canbe considered as a signal of length N/4, that is f′_(T′) ₂ ={f′_(2,0),f′_(2,2), . . . , f′_(2,N/2−2)}. Similarly , since f′₂ _(m) ^(,t) =0 ifg.c.d. (greatest common divisor) (p,t)≠2^(m), consider f′_(T′) ₂ _(m) asthe signal of length N/2^(m+1), and the N/2-point DFT in equation (15)becomes the 2^(r−m−1)-point DFT.

Denoted by Γ′ (f), the totality of splitting-signalsΓ′(f)={f′ _(T′) ₁ , f′ _(T′) ₂ , f′ _(T′) ₄ , . . . , f′ _(T′) ₂ _(r−1), f′ _(T′) ₀ }describes uniquely the signal f in the paired representation. This setof splitting signals is the image of the paired transformation over thesignal f.

As an example, FIG. 3 shows (a) a signal of length 256 and (b) itspaired transform. Eight splitting signals are separated by the verticallines and can be seen in part (b). In Example 1, the calculation andsplitting of a paired transform of a signal is shown.

1-D Linear Convolution Realization

Consider the cyclic convolution of two sequences f_(k) and y_(K) in theframework of the linear filtration. Let Y be a 1-D linear filter withthe impulse characteristics y_(m), and let {circumflex over (f)} be thefiltered signal which is described by the discrete linear cyclicconvolution

$\begin{matrix}\begin{matrix}{{{\hat{f}}_{n} = {\sum\limits_{k = 0}^{N - 1}{y_{n - k}f_{k}}}},} & {{n = {0:( {N - 1} )}},}\end{matrix} & (18)\end{matrix}$where the indices modulo N are taken. The linear filtration in thefrequency domain is described by{circumflex over (F)} _(p) =Y _(p) F _(p) , p=0:(N−1),  (19)where the frequency characteristics Y_(p), and {circumflex over (F)}_(p)and F_(p) are the discrete Fourier transforms of y_(n), and the signals{circumflex over (f)}_(n) and f_(n), respectively.

In terms of the paired representation, the spectral components{circumflex over (F)}_(p) and F_(p) are described by the correspondingsplitting-signals of length N/2{circumflex over (f)}_(T′) _(p) ={{circumflex over (f)}′_(p,0),{circumflex over (f)}′_(p,1), . . . , {circumflex over (f)}′_(p,N/2−1)}{circumflex over (f)}_(T′) _(p) ={f′_(p,0), f′_(p,1), . . . ,f′_(p,N/2−1)}  (20) and (21)and their values are defined by the paired functions as{circumflex over (f)}′ _(p,t)=χ′_(p,t) ∘{circumflex over (f)}, f′_(p,t)=χ′_(p,t) ∘f, t=0:(N/2−1)  (22)

Owing to the construction of the partition σ′ in (17), the pairedtransformation χ′_(N) is induced by the functions χ′_(p,t) for p=2^(n)(n=0:(r−1)) and p=0. Therefore, it is sufficient to analyze the equality(18) only at these r+1 points. The response function Y_(p) of the linearfilter in the polar form (|Y_(p)|, φ_(p)) is considered, where |Y_(p)|is the absolute value and φ_(p)ε[−π, π] is the phase of Y_(p). The phaseφ_(p) by the angle composed by the line drawn from the original andintersected the unit circle at the nearest point of the partition of thecircle divided to N equal parts is approximated. Take the value ofφ_(p)=φ_(p)/(2π/N) and denote by φ _(p)=[φ_(p)] the integer part of thenumber φ_(p). The value of, φ _(p) is a value within the interval [−N/2,N/2].

φ_(p) is assumed to be an integer. This means that the values φ_(p) ofthe filter phases are considered to be quantized in the interval [−N/2,N/2]. Consequently, one can represent the response function of thefilter in the following formY _(p) =|Y _(p) |W ^(φp) , p=2^(n) , n=0:(r−1)  (23)Apparently, for large values N (for instance, such as N≧1024), it can beassumed that the normed phase φ _(p)=[φ_(p)=(2π/N)]≈φ_(p)=(2π/N). Thatis why, for large N, φ _(p)=φ_(p)/(2π/N) (See FIG. 4) may be used.Furthermore for practical applications, when the filter frequencycharacteristics Y_(p) is computed via the N-point DFT, the phases ofY_(p) are exact multiples of 2π/N. The number 2π/N is the common periodof all basis exponential functions W^(np) of the DFT (15). Therefore,the equality in (23) holds for all frequencies p under consideration.

If φ _(p)<0, the right part of (19) can be written as

$\begin{matrix}{{Y_{p}F_{p}} = {{\sum\limits_{t = 0}^{{N/2} - 1}{( {{Y_{p}}f_{p,t}^{\prime}} )W^{t + {\overset{\_}{\varphi}p}}}} = {\sum\limits_{t = 0}^{{N/2} - 1}{( {{Y_{p}}f_{p,{t - {\overset{\_}{\varphi}p}}}^{\prime}} )W^{t}}}}} & (24)\end{matrix}$where f′_(p,−t)=−f′_(p,N/2−t) for t=0:(N/2−1).

If φ _(p)>0, the right part of (19) can be written as

$\begin{matrix}{{Y_{p}F_{p}} = {{\sum\limits_{t = 0}^{{N/2} - 1}{( {{Y_{p}}f_{p,t}^{\prime}} )W^{t - {\overset{\_}{\varphi}p}}}} = {\sum\limits_{t = 0}^{{N/2} - 1}{( {{Y_{p}}f_{p,{t + {\overset{\_}{\varphi}p}}}^{\prime}} )W^{t}}}}} & (25)\end{matrix}$where f′_(p,N/2+t)=−f′_(p,t) for t=0:(N/2−1).

Comparing (24) and (25) with (19),

${\hat{F}}_{p} = {{\sum\limits_{t = 0}^{{N/2} - 1}{{\hat{f}}_{p,t}^{\prime}W^{t}}} = {\sum\limits_{t = 0}^{{N/2} - 1}{( {{Y_{p}}f_{p,{t \mp {\overset{\_}{\varphi}p}}}^{\prime}} )W^{t}}}}$

One can see that the components of the splitting-signal {circumflex over(f)}_(T′) _(p) for the filtered image are defined as

$\begin{matrix}{{\hat{f}}_{p,t}^{\prime} = \{ \begin{matrix}{{{Y_{p}}f_{p,{t + {{{sgn}{({\overset{\_}{\varphi}}_{p})}}{\overset{\_}{\varphi}}_{p}}}}^{\prime}},} & {{{if}\mspace{14mu} 0} \leq {t + {{{sgn}( {\overset{\_}{\varphi}}_{p} )}{\overset{\_}{\varphi}}_{p}}} < {N/2}} \\{{{- {Y_{p}}}f_{p,{{({t + {{{sgn}{({\overset{\_}{\varphi}}_{p})}}{\overset{\_}{\varphi}}_{p}}})}\mspace{11mu}{mod}\mspace{11mu}{N/2}}}^{\prime}},} & {otherwise}\end{matrix} } & (26)\end{matrix}$where t=0:(N/2−1), and sgn(x) is the function defined as sgn(x)=−1 ofx<0, and sgn(x)=1 if x≧0.

To calculate the splitting-signal (21) of the filtered signal, one canmultiply the splitting-signal (20) by the scalar |Y_(p)| and shiftcyclicly each of its components f′_(p,t) by φ _(p) elements. Theshifting is to the right, if φ _(p)<0, or to the left, if φ _(p)>0.

Let {circumflex over (F)}_(p) and F_(p) be respectively thesplitting-signals {circumflex over (f)}_(T′) _(p) and f_(T′) _(p)written in the vector forms. Then, in the paired representation, one canwrite equation (19) in the following matrix form{circumflex over (F)} _(p) =|Y _(p) |T ^(F) ^(p) F _(p)  (27)where T is the operator of the cyclic shift of the vector of dimensionN/2 by one element to the right (if φ _(p)<0)TF _(p)=(−{circumflex over (f)}′_(p,N/2−1) , {circumflex over (f)}′_(p,0) , {circumflex over (f)}′ _(p,1) , . . . , {circumflex over (f)}′_(p,N/2−2))′,  (28)and T ^(φ) ^(p) is its φ _(p)-th power.

The matrix of the operator T ^(φ) ^(p) has the form

${T^{{\overset{\_}{\varphi}}_{p}}} = {\begin{matrix}{ \begin{matrix}0 & . & . & 0 & {- 1} & 0 & . & . & . & 0 \\0 & 0 & . & . & 0 & {- 1} & 0 & . & . & 0 \\0 & 0 & 0 & . & . & 0 & {- 1} & 0 & . & 0 \\. & . & . & . & . & . & . & . & . & . \\. & . & . & . & . & . & . & . & 0 & {- 1}\end{matrix} \}{\overset{\_}{\varphi}}_{p}} \\\begin{matrix}1 & 0 & . & . & . & . & . & . & . & . \\0 & 1 & 0 & . & . & . & . & . & . & . \\. & . & . & . & . & . & . & . & . & . \\0 & . & . & 1 & 0 & . & . & . & . & .\end{matrix}\end{matrix}}$if φ _(p)<0, and [T ^(φ) ^(p) ] is described by the matrix transposed to[T^(− φ) ^(p) ] if φ _(p)>0. Note, that [T ^(φ) ^(p) ]=−[T^(N/2− φ) ^(p)].

For the considered points p=2^(n), (n=0: (r−1)), and p=0, the componentsf′_(p,t) and {circumflex over (f)}_(p,t)≡0 in vectors (20) and (21), ifg.c.d.(p,t)≠p. Therefore, the values φ _(p) in (27) can be considered tobe equal to φ _(p)/p, and the size of the matrix [T ^(φ) ^(p) ] to beequal to N/(2p)×N/(2p).

To obtain the estimation {circumflex over (f)} of the initial signal,the inverse discrete paired transform must be fulfilled. The diagram ofrealization the linear convolution by the paired transforms is given inFIG. 1 b.

Linear Convolution in Matrix Form

The simple matrix representation for the convoluted sequence f_(n) withy_(n) is given here.

denotes the operation of the direct product of matrices, and by{circumflex over (f)} and f the vectors composed of signals {circumflexover (f)}_(n) and f_(n), respectively.

Let Y_(p), p=0: (N−1), be the Fourier transform (or, the frequencycharacteristic of the linear filter). Assume that for the phase ofY_(p), the following holds φ _(p)=φ_(p)=(2π/N), for p=0, 1, 2, 4, 8, . .. , N/2. Then, {circumflex over (f)}_(n) in can be written in thefollowing matrix form

$\begin{matrix}{\hat{f} = {{( \chi_{N}^{\prime} )^{- 1}}( {\underset{n = 0}{\overset{r}{\otimes}}{( {{Y_{2^{n}}} \otimes {I_{2^{r - n - 1}}}} )\lbrack T^{{\overset{\_}{\varphi}}_{2^{n}}} \rbrack}} ){\chi_{N}^{\prime}}{f.}}} & (29)\end{matrix}$I_(M) is the unit matrix M×M.

For simplification purposes, consider |Y_(2r)| to be equal to |Y₀|, φ ₂_(r) = φ ₀≡0, and I⁻¹=1:

Estimation of Arithmetical Operations

Let N=2^(r), r>>1, then for the realization of the linear convolution ofthe real sequence f_(n) with y_(n) by means of the paired transform itis sufficiently to fulfillM′ _(N) =N, and C′ _(N)=4N−4  (30)real operations of multiplications and additions, respectively. In fact,from the matrix representation (29) it follows that for the describedfiltering it is necessary to fulfill N real multiplications by r+1modules of the frequency characteristics Y_(p). The number of additionsrequired for the direct or the inverse paired transformation equals2N−2.

For a complex signal {circumflex over (f)}, the number of realoperations of multiplications equals 2N, since two real multiplicationsare needed to multiply a complex number by a real one (modulus). Thenumber of additions in this case equals 2(4N−4). The number of thearithmetical operations is exceeded in real case by two times.

Thus, the described method requires O(N) arithmetical operations. Itshould be noted for the comparison, that the traditional method of theDFT uses O(rN) arithmetical operations of multiplications and additions.Indeed, if considering that the number of the complex multiplicationsfor the N-point DFT equals 2^(r−1)(r−3)+2, then the number of complexmultiplications required for filtering the signal via the DFT methodequals2|2^(r−1)(r−3)+2|+2^(r)=2^(r)(r−2)+4.In a case where one complex multiplication is executed via three realmultiplications and three real additions, then the numbers of realmultiplications and additions required to calculate the convolution bythe FFT method can be estimated as follows:M _(N) ¹=3|2^(r)(r−2)+4|=3×2^(r)(r−2)+12A _(N)¹=3|2(2^(r−1)(r−3)+2)+r2^(r)+2^(r)|+4|2^(r)−2|=3×2^(r)(r+3)+4.  (31)Taking into account the property of complex conjugation of spectralcomponents for the real data, {circumflex over (F)}*_(p)={circumflexover (F)}*_(N−p) for p=1: (N/2−1), the following estimates are obtainedM _(N) ¹=3/2N(r−2)+6, A _(N) ¹=3/2N(r+3)+2.  (32)

The proposed method uses N multiplications, i.e. it requires μ_(N) timesless multiplications than the fast DFT method uses, where

$\mu_{N} = {\frac{{3/2}{N( {r - 2} )}}{N} = {{3/2}{( {r - 2} ).}}}$

For example, μ₂₅₆=9, μ₁₀₂₄=12, and μ₄₀₉₆=15. If compare the operationsof additions, then the advantage of the proposed algorithm is estimatedby coefficient

$\alpha_{N} = {\frac{{3{N/2}( {r + 3} )} + 4}{{4N} - 4} \geq {{3/8}{( {r + 3} ).}}}$

For example, α₂₅₆≈4, α₁₀₂₄≈5, and α₄₀₉₆≈5.6.

Only (r+1) values of the Fourier transform Y_(p) (or, frequencycharacteristics of the filter) are used in the proposed method.Additional operations reduction can be achieved by using a non completeDFT algorithm for computing required Y_(p).

Calculation of Inverse Paired Transform

Unitary Inverse Paired Transform

Consider the transform that is inverse to the unitary paired transformχ′_(N). The N-point discrete inverse unitary transform (χ′)_(N) ⁻¹ isdefined by the matrix Y_(N) ¹ being inverse to the matrix X_(N) ¹ of theN-point discrete inverse unitary transform, i.e.Y¹ _(N)X¹ _(N)=I_(N)where I_(N) is the unit matrix N×N.

The matrix of the inverse paired transform can be easily found if weconsider the paired transform with normalized basis functions. In thiscase, the inverse paired transform is defined by the matrix Y_(N) ¹.that is transposed to the matrix X_(N) ¹.

Inverse 8-Point DPT

Consider the matrix of the 8-point paired transform χ′₈

$X_{8}^{1} = \begin{bmatrix}1 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 \\0 & 1 & 0 & 0 & 0 & {- 1} & 0 & 0 \\0 & 0 & 1 & 0 & 0 & 0 & {- 1} & 0 \\0 & 0 & 0 & 1 & 0 & 0 & 0 & {- 1} \\1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} & 0 \\0 & 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} \\1 & {- 1} & 1 & {- 1} & 1 & {- 1} & 1 & {- 1} \\1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\end{bmatrix}$

The paired transform of the signal f={0, 1, 2, 3, 4, 5, 6, 7} is thesignal g={−4, −4, −4, −4, −4, −4, −4, 28}.

Indeed

${\begin{bmatrix}1 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 \\0 & 1 & 0 & 0 & 0 & {- 1} & 0 & 0 \\0 & 0 & 1 & 0 & 0 & 0 & {- 1} & 0 \\0 & 0 & 0 & 1 & 0 & 0 & 0 & {- 1} \\1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} & 0 \\0 & 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} \\1 & {- 1} & 1 & {- 1} & 1 & {- 1} & 1 & {- 1} \\1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\end{bmatrix}\begin{bmatrix}0 \\1 \\2 \\3 \\4 \\5 \\6 \\7\end{bmatrix}} = \begin{bmatrix}{- 4} \\{- 4} \\{- 4} \\{- 4} \\{- 4} \\{- 4} \\{- 4} \\28\end{bmatrix}$

The inverse 8-point paired transform has the matrix

$Y_{8}^{1} = {\frac{1}{8}\begin{bmatrix}1 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 \\0 & 1 & 0 & 0 & 0 & {- 1} & 0 & 0 \\0 & 0 & 1 & 0 & 0 & 0 & {- 1} & 0 \\0 & 0 & 0 & 1 & 0 & 0 & 0 & {- 1} \\1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} & 0 \\0 & 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} \\1 & {- 1} & 1 & {- 1} & 1 & {- 1} & 1 & {- 1} \\1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\end{bmatrix}}$

The inverse transform of the obtained g={−4, −4, −4, −4, −4, −4, −4, 28}results in the initial singal f={0, 1, 2, 3, 4, 5, 6, 7}. Indeed,

${{\frac{1}{8}\begin{bmatrix}1 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 \\0 & 1 & 0 & 0 & 0 & {- 1} & 0 & 0 \\0 & 0 & 1 & 0 & 0 & 0 & {- 1} & 0 \\0 & 0 & 0 & 1 & 0 & 0 & 0 & {- 1} \\1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} & 0 \\0 & 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} \\1 & {- 1} & 1 & {- 1} & 1 & {- 1} & 1 & {- 1} \\1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\end{bmatrix}}\begin{bmatrix}{- 4} \\{- 4} \\{- 4} \\{- 4} \\{- 4} \\{- 4} \\{- 4} \\28\end{bmatrix}} = \begin{bmatrix}0 \\1 \\2 \\3 \\4 \\5 \\6 \\7\end{bmatrix}$

Inverse 16-Point DPT

Consider the matrix of the 16-point paired transform χ′₁₆

$X_{16}^{1} = \begin{bmatrix}1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 1} & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 1} & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 1} \\1 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 & 1 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 \\0 & 1 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 & 1 & 0 & 0 & 0 & {- 1} & 0 & 0 \\0 & 0 & 1 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 & 1 & 0 & 0 & 0 & {- 1} & 0 \\0 & 0 & 0 & 1 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 & 1 & 0 & 0 & 0 & {- 1} \\1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} & 0 \\0 & 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} \\1 & {- 1} & 1 & {- 1} & 1 & {- 1} & 1 & {- 1} & 1 & {- 1} & {- 1} & {- 1} & {- 1} & {- 1} & {- 1} & {- 1} \\1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & {- 1} & 1 & {- 1} & 1 & {- 1} & 1 & {- 1}\end{bmatrix}$

The 16-point paired transform of the signal f={2, 3, 2, 5, 2, 3, 5, 7,1, 2, 4, 8, 7, 4, 2, 1} results in the output g=X₁₆ ¹; f′={1, 1, −2, −3,−5, −1, 3, 6, −6, −2, −1, 5, −1, −9, −8, 58}′.

The inverse 16-point paired transform has the matrix

$X_{16}^{1} = {\frac{1}{16}\begin{bmatrix}8 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 8} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 8 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 8} & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 8 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 8} & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 8 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 8} & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 8 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 8} & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 8 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 8} & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 8 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 8} & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 8 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 8} \\4 & 0 & 0 & 0 & {- 4} & 0 & 0 & 0 & 4 & 0 & 0 & 0 & {- 4} & 0 & 0 & 0 \\0 & 4 & 0 & 0 & 0 & {- 4} & 0 & 0 & 0 & 4 & 0 & 0 & 0 & {- 4} & 0 & 0 \\0 & 0 & 4 & 0 & 0 & 0 & {- 4} & 0 & 0 & 0 & 4 & 0 & 0 & 0 & {- 4} & 0 \\0 & 0 & 0 & 4 & 0 & 0 & 0 & {- 4} & 0 & 0 & 0 & 4 & 0 & 0 & 0 & {- 4} \\2 & 0 & {- 2} & 0 & 2 & 0 & {- 2} & 0 & 2 & 0 & {- 2} & 0 & 2 & 0 & {- 2} & 0 \\0 & 2 & 0 & {- 2} & 0 & 2 & 0 & {- 2} & 0 & 2 & 0 & {- 2} & 0 & 2 & 0 & {- 2} \\1 & {- 1} & 1 & {- 1} & 1 & {- 1} & 1 & {- 1} & 1 & {- 1} & {- 1} & {- 1} & {- 1} & {- 1} & {- 1} & {- 1} \\1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & {- 1} & 1 & {- 1} & 1 & {- 1} & 1 & {- 1}\end{bmatrix}}$

One can check that the inverse transform of the signal g={1, 1, −2, −3,−5, −1, 3, 6, −6, −2, −1, 5, −1, −9, −8, 58} the initial signal f={2, 3,2, 5, 2, 3, 5, 7, 1, 2, 4, 8, 7, 4, 2, 1}, i.e.X₁₆ ¹g′=f′.

The methods described previously may be implemented as a program that isrun on a computer, or by a chip or microcontroller. These methods mayalso be run on other programs such as MatLab (commercially availablefrom The MathWorks, Inc., Natick, Mass.). Examples 4 and 5 show examplesof MatLab programs which would implement these methods.

The terms a or an, as used herein, are defined as one or more than one.The term plurality, as used herein, is defined as two or more than two.The term means, as used herein, is defined as hardware, firmware and/orsoftware for achieving a result. The term program or phrase computerprogram, as used herein, is defined as a sequence of instructionsdesigned for execution on a computer system. A program, or computerprogram, may include a subroutine, a function, a procedure, an objectmethod, an object implementation, an executable application, an applet,a servelet, a source code, an object code, a shared library/dynamic loadlibrary and/or other sequence of instructions designed for execution ona computer system. The phrase any integer derivable therein, as usedherein, is defined as an integer between the corresponding numbersrecited in the specification.

EXAMPLES

Specific embodiments of the invention will now be further described bythe following, nonlimiting examples which will serve to illustrate insome detail various features. The following examples are included tofacilitate an understanding of ways in which the invention may bepracticed. It should be appreciated that the examples which followrepresent embodiments discovered to function well in the practice of theinvention, and thus can be considered to constitute preferred modes forthe practice of the invention. However, it should be appreciated thatmany changes can be made in the exemplary embodiments which aredisclosed while still obtaining like or similar result without departingfrom the spirit and scope of the invention. Accordingly, the examplesshould not be construed as limiting the scope of the invention.

Example 1

Let N=16, and let the signal be f={2, 3, 2, 5, 2, 3, 5, 7, 1, 2, 4, 8,7, 4, 2, 1}. Then, for points p=1, 2, 4, 8, and 0, the correspondingfive splitting-signals of Γ′ (f) are defined as follows

$\begin{matrix}\begin{matrix}{f_{T_{1}} = \{ {{f_{0} - f_{8}},{f_{1} - f_{9}},{f_{2} - f_{10}},{f_{3} - f_{11}}} } \\ {{f_{4} - f_{12}},{f_{5} - f_{13}},{f_{6} - f_{14}},{f_{7} - f_{15}}} \} \\{= \{ {1,1,{- 2},{- 3},{- 5},{- 1},3,6} \}}\end{matrix} \\\begin{matrix}{f_{T_{2}} = \{ {{f_{0} - f_{4} + f_{8} - f_{12}},{f_{1} - f_{5} + f_{9} - f_{13}}} } \\ {{f_{2} - f_{6} + f_{10} - f_{14}},{f_{3} - f_{7} + f_{11} - f_{15}}} \} \\{= \{ {{- 6},{- 2},{- 1},5} \}}\end{matrix} \\\begin{matrix}{f_{T_{4}} = \{ {f_{0} - f_{2} + f_{4} - f_{6} + f_{8} - f_{10} + f_{12} - f_{14}} } \\ {f_{1} - f_{3} + f_{5} - f_{7} + f_{9} - f_{11} + f_{13} - f_{15}} \} \\{= \{ {{- 1},{- 9}} \}}\end{matrix} \\{f_{T_{5}} = {{f_{0} - f_{1} + \ldots + f_{14} - f_{15}} = {- 8}}} \\{f_{T_{a}} = {{f_{0} - f_{1} + \ldots + f_{14} - f_{15}} = 58}}\end{matrix}$

And they compose the 16-point paired transform

$\begin{matrix}{{X_{16}^{\prime} \diamond \; f} = \{ \underset{︸}{1,1,{- 2},{- 3},{- 5},{- 1},3,6} } \\{\underset{︸}{{- 6},{- 2},{- 1},5}} \\{\underset{︸}{{- 1},{- 9}}} \\{\underset{︸}{- 8}} \\ \underset{︸}{58} \}\end{matrix}$

Example 2

A 16-point convolution is shown here, in accordance with an embodimentof the present disclosure.

Let N=16 and let the impulse characteristic of the filter bey={0, . . . , 0, 1/11, 2/11, 5/11, 2/11, 1/11, 0, . . . 0}

with value y₈=5/11 in the middle. The 16-point DFT of y(n) at point p=1isY₁=−0.8491−j0.3517=|0.9191|exp(−j2.7489).

So, |Y₁|=0.9191 and φ₁=−2.7489 (in radians). Since, 2π/16=0.3927, thefollowing normalized phase may be obtainedφ ₁=φ₁/0.3927=−7.

The frequency characteristic at point 1 can be written asY ₁=0.9191W ⁷ ·W=W ₁₆=exp(−j2π/16).Letf_(n)={2, 3, 2, 5, 2, 3, 5, 7, 1, 2, 4, 8, 7, 4, 2, 1}be the 16-point signal to be processed by the filter. Thesplitting-signal that corresponds to the point p=1

$\begin{matrix}{f_{T_{1}}^{\prime} = \{ {f_{1,0}^{\prime},f_{1,1}^{\prime},f_{1,2}^{\prime},f_{1,3}^{\prime},f_{1,4}^{\prime},f_{1,5}^{\prime},f_{1,6}^{\prime},f_{1,7}^{\prime}} \}} \\{= \{ {1,1,{- 2},{- 3},{- 5},{- 1},3,6} \}}\end{matrix}$should be shifted by, according to the operator T⁷ that has the 8×8matrix,

$T^{7} = {\begin{bmatrix}\; & {- 1} & \; & \; & \; & \; & \; & \; \\\; & \; & {- 1} & \; & \; & \; & \; & \; \\\; & \; & \; & {- 1} & \; & \; & \; & \; \\\; & \; & \; & \; & {- 1} & \; & \; & \; \\\; & \; & \; & \; & \; & {- 1} & \; & \; \\\; & \; & \; & \; & \; & \; & {- 1} & \; \\\; & \mspace{11mu} & \; & \; & \; & \; & \; & {- 1} \\1 & \; & \; & \; & \; & \; & \; & \;\end{bmatrix} = {- ( T^{1} )^{T}}}$seven points to the right

$\begin{matrix} f_{T_{1}}^{\prime}arrow\{ {{- f_{1,1}^{\prime}},{- f_{1,2}^{\prime}},{- f_{1,3}^{\prime}},{- f_{1,4}^{\prime}},{- f_{1,5}^{\prime}},{- f_{1,6}^{\prime}},{- f_{1,7}^{\prime}},f_{1,0}^{\prime}} \}  \\{= \{ {{- 1},{+ 2},{+ 3},{+ 5},{+ 1},{- 3},{- 6},1} \}}\end{matrix}$Therefore, the splitting-signal of the filtered signal is represented as{circumflex over (f)}′ _(T) ₁ =|Y ₁|{−1, 2, 3, 5, 1, −3, −6, 1}.The DFT of the degraded signal at point p=1 is calculated as{circumflex over (F)} ₁ =|Y ₁||(−1)W ₁₆ ⁰+(+2)W ₁₆ ¹+(+3)W ₁₆ ²+(+5)W ₁₆³+(+1)W ₁₆ ⁴+(−3)W ₁₆ ⁵+(−6)W ₁₆ ⁶+(+1)W ₁₆ ⁷|Indeed, the DFT of the signal f at point p=1 isF ₁=|(+1)W ₁₆ ⁰+(+1)W ₁₆ ¹+(−2)W ₁₆ ²+(−3)W ₁₆ ³+(−5)W ₁₆ ⁴+(−1)W ₁₆⁵+(+3)W ₁₆ ⁶+(6)W ₁₆ ⁷|Multiplying F₁ by Y₁=|Y₁|W₁₆ ⁷, one can obtain

$\begin{matrix}{{Y_{1}F_{1}} = {{Y_{1}}{{{( {+ 1} )W_{16}^{7}} + {( {+ 1} )W_{16}^{8}} + {( {- 2} )W_{16}^{9}} + {( {- 3} )W_{16}^{10}} +}}}} \\{{\mspace{31mu}{{( {- 5} )W_{16}^{11}} + {( {- 1} )W_{16}^{12}} + {( {+ 3} )W_{16}^{13}} + {(6)W_{16}^{14}}}} =} \\{{\hat{F}}_{1} = {{Y_{1}}{{{( {+ 1} )W_{16}^{7}} + {( {- 1} )W_{16}^{0}} + {( {+ 2} )W_{16}^{1}} + {( {+ 3} )W_{16}^{2}} +}}}} \\{\mspace{31mu}{{( {+ 5} )W_{16}^{3}} + {( {+ 1} )W_{16}^{4}} + {( {- 3} )W_{16}^{5}} + {( {- 6} )W_{16}^{6}}}}\end{matrix}$Thus, the DFT of the degraded signal at points p=2k+1, k=1:(N/2−1), iscalculated as{circumflex over (F)} ₂₆₊₁ =|Y ₁||(−1)+(2W ₁₆ ¹)W ₈ ^(k)+(3W ₁₆ ²)W ₈^(2k)+(5W ₁₆ ³)W ₈ ^(3k)+(1W ₁₆ ⁴)W ₈ ^(4k)+(−3W ₁₆ ⁵)W ₈ ^(5k)+(−6W ₁₆⁶)W ₈ ^(6k)+(1W ₁₆ ⁷)W ₈ ^(7k)|.One can process similarly the splitting-signal, if the phase φ ₁ takesanother values. For instance, if the phase φ ₁=−3, then the sequencef′_(T1) will be shifted (to the right) by operator T³ that has thematrix

$T^{3} = \begin{bmatrix}\; & \; & \; & \; & \; & {- 1} & \; & \; \\\; & \; & \; & \; & \; & \; & {- 1} & \; \\\; & \; & \; & \; & \; & \; & \; & {- 1} \\1 & \; & \; & \; & \; & \; & \; & \; \\\; & 1 & \; & \; & \; & \; & \; & \; \\\; & \; & 1 & \; & \; & \; & \; & \; \\\; & \; & \; & 1 & \; & \; & \; & \; \\\; & \; & \; & \; & 1 & \; & \; & \;\end{bmatrix}$If the phase φ ₁=+3, then the sequence f′_(T1) will be shifted (to theleft) by operator T³ that has the matrix

$T^{- 3} = \begin{bmatrix}\; & \; & \; & {1\;} & \; & \; & \; & \; \\\; & \; & \; & \; & {\; 1} & \; & \mspace{11mu} & \; \\\; & \; & \; & \; & \; & {\; 1} & \; & \; \\\; & \; & \; & \; & \; & \; & {1\;} & \; \\\; & \; & \; & \; & \; & \; & \; & {1\;} \\{\;{- 1}} & \; & \; & \; & \; & \; & \; & \; \\\; & {\;{- 1}} & \; & \; & \; & \; & \; & \; \\\; & \; & {{- 1}\;} & \; & \; & \; & \; & \;\end{bmatrix}$Note, that T⁻³=(T³)^(T)=−T⁵.

Example 3

A 256-point convolution is shown here, in accordance with an embodimentof the present disclosure.

Consider the signal f(t) of length N=256 shown in FIG. 3( a) and thesmoothing filter with the impulse characteristicy={0, . . . , 0, 1, 2, 5, 2, 1, 0, . . . 0}/11shifted to the center N/2+1=129.

At nine points p=0 and p=2^(n), n=0:7, the response function of thesmoothing filter is represented by the values from the following table

p 0 1 2 4 8 16 32 64 128 |Y_(p)| 1 0.9997 0.9987 0.9948 0.9792 0.01010.7117 0.2727 0.2727 φ_(p) 0 −3.117 0.0491 0.0982 0.1983 0.3927 0.78541.5708 3.1416 φ _(p) 0 −127 2 4 8 16 32 64 128The paired transform χ′[f] of the signal as well as its splitting into(log₂ N+1 =9) short signalsχ′|f|={f′_(T) ₁ , f′_(T) ₂ , f′_(T) ₄ , f′_(T) ₈ , f′_(T) ₁₆ , f′_(T) ₃₂, f′_(T) ₆₄ , f′_(T) ₁₂₈ , f′_(T) ₀ }are shown in FIG. 3( b).

According to this table, the short signals are processed in thefollowing way: Step 1: (p=1) The first signal of length 128 is shiftedcyclicly by 127 elements to the right (or, 1 element to the left), thesigns of the first 127 components are changed, and all components aremultiplied by |Y₁|=0.9997.f′ _(T) ₁ ={f′ _(1,0) , f′ _(1,1) , f′ _(1,2) , . . . , f′ _(1,125) , f′_(1,126) , f′ _(1,127) }→{f′ _(1,1) , f′ _(1,2) , . . . , f′ _(1,125) ,f′ _(1,126) , f′ _(1,127) , f′ _(1,0) }→{−f′ _(1,1) , −f′ _(1,2) , . . ., −f′ _(1,125) , −f′ _(1,126) , −f′ _(1,127) , f′ _(1,0)}×0.9997.Step 2: (p=2) The second signal of length 64 is shifted cyclicly by 1(because φ ₂/2=1) element to the left, the sign of the last component ischanged, and all components are multiplied by |Y₂|=0.9987:f′ _(T) ₂ ={f′ _(2,0) , f′ _(2,2) , f′ _(2,4) , . . . , f′ _(2,122) , f′_(2,124) , f′ _(2,126) }→{f′ _(2,2) , f′ _(2,4) , . . . , f′ _(2,122) ,f′ _(2,124) , f′ _(2,126) , f′ _(2,0) }→{f′ _(2,2) , f′ _(2,4) , . . . ,f′ _(2,122) , f′ _(2,124) , f′ _(2,126) , −f′ _(2,0)}×0.9987.Step 3: (p=4) The third signal of length 32 is shifted cyclicly by 1(because φ ₄/4=1) element to the left, the sign of the last component ischanged, and all components are multiplied by |Y₄|0.9948:f′ _(T) ₄ ={f′ _(4,0) , f′ _(4,4) , f′ _(4,8) , . . . , f′ _(4,116) , f′_(4,120) , f′ _(4,124) }→{f′ _(4,4) , f′ _(4,8) , . . . , f′ _(4,116) ,f′ _(4,120) , f′ _(4,124) , f′ _(4,0) }→{f′ _(4,4) , f′ _(4,8) , . . . ,f′ _(4,116) , f′ _(4,120) , f′ _(4,124) , −f′ _(4,0)}×0.9948.Steps 4-7: (p=8, 16, 32, 64) The next four signals of length 16, 8, 4,and 2 are shifted cyclicly by 1 (because φ _(p)/p=1) element to theleft, the sign of the last components are changed, and all componentsare multiplied by |Y_(p)| as follows:f′_(T) ₈ ={f′ _(8,0) , f′ _(8,8) , f′ _(8,16) , . . . , f′ _(8,104) , f′_(8,112) , f′ _(8,120) }→{f′ _(8,8) , f′ _(8,16) , . . . , f′ _(8,104) ,f′ _(8,112) , f′ _(8,120) , −f′ _(8,0)}×0.9792,f′ _(T) ₁₆ ={f′ _(16,0) , f′ _(16,16) , f′ _(16,32) , . . . , f′_(16,94) , f′ _(16,112) }→{f′ _(16,16) , f′ _(16,32) , . . . , f′_(16,94) , f′ _(16,112) , −f′ _(16,0)}×0.9191,f′ _(T) ₃₂ ={f′ _(32,0) , f′ _(32,32) , f′ _(32,64) , f′ _(32,96) ,}→{f′ _(32,32) , f′ _(32,64) , f′ _(32,96) , −f′ _(32,0)}×0.7117,f′ _(T) ₆₄ ={f′ _(64,0) , f′ _(64,64) }→{f′ _(64,64) , −f′_(64,0)}×0.2727.Step 8: (p=128) The one-point signal f′_(T) ₁₂₈ is processed asf′ _(T) ₁₂₈ ={f′ _(128,0) }→{−f′ _(128,0)}×0.2727.Step 9: (p=0) The last one-point signal f′_(T) ₀ is processed asf′ _(T) ₀ ={f′ _(0,0) }→{f′ _(0,0)}×1.The results of the above processing may be denoted byχ′|g|={g′ _(T) ₁ , g′ _(T) ₂ , g′ _(T) ₄ , g′ _(T) ₈ , g′ _(T) ₁₆ , g′_(T) ₃₂ , g _(T) ₆₄ , g _(T) ₁₂₈ , g _(T) ₀ }which is the paired transform of the filtered signal. FIG. 5( d) showsthe transform χ′[g] composed by the processed 9 signals. The inversepaired transform of χ′[g], or the filtered signal is shown in FIG. 5(c). The spectral characteristics of absolute value and phase of thefilter, as well as the values that have been used (marked by circles)are shown in FIG. 6.

Example 4

Shown here is an embodiment of a MatLab implementation of convolutionmethods, in accordance with embodiments of this disclosure.

% Call: dcc_realization.m % Dr. Art M. Crigoryan, EX UTSA 03/09-23/2002% % 1. This is the demo program to illustrate the realization % of thelinear cyclic convolution (or linear filtration) % %    g(n)=f(n)*h|n|−sum _(k−1){circumflex over ( )}256 f(n−k)h|k| %       n, k=1,2,3, ..., 256, % % of signal f|n| with the sequence h|n|being the impulse % characteristic of the smoothing filter, which hasthe % form of the triangle (1 2 5 2 1) (normalized). % The program isbased on the fast algorithms of the direct % and inverse pairedtransforms, (fpt_1d{} and fpt1_1d{}). % The Fourier transform H(w) ofthe impulse characteristics % is computed by the fast Fourier transform{fft()} % %     H[w}=sum_(n=0){circumflex over ( )}126f(n)exp(=j2pi/256*nw) %        w=0,1,2, ..., 255. % % The code user only9 values of H{w}, namely at points % w=1,2,4,8,16,32,64,128, and 0. % %2. For comparison with the FFT method (that uses all 256 % values ofH[w]) the result of applying the FFT for the % convolution realizationin also illustrated. % % 3. This program used the 1-D signal of length256 written % in the file ‘image_signal.sig’ % % 4. Used codes: %  fpt_ld{} -- fast paired transform [Artycom] %   fpti_ld{} -- inversefast paired transform [Artyom] %1. %A. Reading the signal from file------------------------ fid=fopen{‘image_signal,sig’ , ‘rb’};f=fread(fid, 256, ‘double’);  f=f′;  fclose(fid): clear fid:N=length(f):  % N=256 % open the figure and display the original signalf_dft1−figure(‘Position’, (10 398 560 420), ...   ‘Tag’, ‘FigD’, ...  ‘Name’,‘ 1-D Signal decomposition v.1, Art M.G., 2002’);subplot(2,2,1); plot(f); text (230,158,‘t’); b_title=title(‘Originalsignal of length 256’); set(h_title, ‘FontName’, ‘Times’, ‘FontSize’,9); axis([0,260,160,185]); %B. Setting the impulse characteristic h(n)for 1-D filter %  The triangle impulse of length 5 shifted to the centerh1=ones(1,5); h1=[1 2 5 2 1]; s°sum(h1); h1=h1/s;  % normalyzing theimpulse n1=length(h1): L=(n1−1)/2; LL=N/2; h=zeros(1,N);h(LL−L:LL+L)=h1; % triangle is shifted to the center H=fft(b,N); % theN-point DFT of the sequence h(n) %C. Compute log(N)+1 magnitudes andphases components of % transfer function that will be used by FPTmethod: r1=log2 (N)+1; H_only=zeroB(l,r1);  % for r1 absolute values ofthe filter P_only=zeros(l,r1);  % for r1 phase values of the filter p=1;p_need(1)=p; for k=2:r1   p_need(k)=p+1;   p=p*2; end;%p_need=[1,2,3,5,9,17,33,65,129]; H_only=abs(H(p_need));P_only=angle(H(p_need)); %D. Normalizing the angles in P_only % Phaserepresentation by dividing the unit circle by N parts pn=2*pi/N;P_only=round(P_only/pn); fprintf(‘ %2.0f’,P_only); fprintf(‘ \n ’); % 0−127 2 4 8 16 32 64 128 fprintf(‘%6.4f ’,H_only); % 1.0000 0.9997 0.99870.9948 0.9792 0.9191 0.7117 0.2727 0.2727 %2. Performing theconvolution - - - - - - - - %A. Compute the paired transform of thesignal f(n) SP=fpt_ld(f);  % the N=point DPT of the signal %B. Computingthe DPT of the convoluted signal g(n) SP_new=zeros(1,N); % for theoutput DPT %a) computing two last components of the DPT of g(n): p2−1;SP_new(N) = SP(N)*H_only(1); % component #N N2−N/2; T=1; ifP_only(k)==N2 T=−1; end; SP_new(N−1)= T* (SP(N−1))*H_only(r1); %component #N−1 %b) computing the rest (N−2) components of the DPT ofg(n): n_1=1; n_2=N2; for k−2:log2(N}   ph1-P_only(k);  GT=zeros(1,N2};   % for split-signals or the DPT of f(n)  GT(1:p2;N2)=SP(n_1:n_2):   GP = zeros(1, N2);  % for split-signals ofthe DPT of g(n)   M2=N2/p2;   if phi < 0     for m=1:p2:−phi       GP(m)= −GT(M2+phi+m);     end;     m1−m;     mm=1;     for m=m1+p2:N2      GP(m) −GT|mm);       mm=mm+1;     end;   elseif phi >= 0     form=1:p2:N2−phi       GP(m)=GT(ph1+m);     end;     m1=m;     mm=1;    for n=m1+p2:p2:N2       GP(m) = −GT(mn);       nm=mm+1;     end;  end;   SP_new(n_1:n_2)=GP(1:p2:N2) *H_only(k);   p2*p2*2;  n_1=n_2+1;  n_2=n_2+M2/2; end; % To print DPTs of the original andfiltered signals fprintf(‘ %6.4f,’,SP); fprintf(‘ \n ’); fprintf(‘%6.4f,’,SP_new); %e) Displaying the results of the linear convolutionsubplot (2,2,2); plot(SP); axis([0,260, −25,25]); S_T=sprintf(‘1-D DPTof the signal f(t)’); h_title=title(S_T);set(h_title,‘Color’,‘k’,‘FontName’,‘Times’,‘FontSize’,9);subplot(2,2,4); plot(SP_new); %plot(fftshift(SP_new)); axis([0,260,−25,25]); % hold on; plot(−SP, ‘r’); S_T=sprintf(‘1-D DPT of thefiltered signal’); h_title=title(S_T);set(h_title,‘Color’,‘k’,‘FontName’,‘Times’,‘FontSize’,10);f_new=fpti_ld(SP_new); subplot(2,2,3); plot(fftshift|f_new));axis([0,260,160,185]); hold on; S_T=sprintf(‘The signal filtered by theDPT’); h_title=title[S_T];set(h_title,‘Color’,‘k’,‘FontName’,‘Times’,‘FontSize’,9); % print -dpsfilterPT_1.ps % For comparison with the fast DFT method of theconvolution F*fft(f,N);    % the N-point DFT of the original signalFC=F * H; fc=ifft(FC); fc_real=real(fftshift(fc)); plot(fc_real,‘r’); %print -dpsc filterPT_1b.ps %3. Draw thefilter - - - - - - - - - - - - - - - - - - - -f_dpt2-figure(‘Position’,[10 398 560 420], ...   ‘Tag’,‘FigD’, ...  ‘Name’,‘ 1-D Signal decomposition v.1. Art M.G., 2002’); %a) Plot thesbs. of the frequency characteristics of the filter h_1=subplot(2,1,1);plot(sbs(H)); axis([−5,135,0,1.1]);    %    h=[1 2 5 2 ]1/11S_T=sprintf|‘|1-D DFT| of the impulse reponse function’);h_title=title(S_T);set[h_title,‘Color’,‘k’,‘FontName’,‘Times’,‘FontSize’,10]; %b) Mark thevalue that were used by the FPT method hold on; for k=1:r1  x1=p_need(k);   h_o=plot(x1,H_only(k),‘o’);   set(h_o,‘Color’,[1 00]);   x2=[x1 x1], y2=[0,H_only(k)];  h_line1=line(x2,y2,‘LineStyle’,‘−’);   set(h_line1,‘Color’,[1 0 0]);end; %e) Plot the phase of the filter h_2=subplot[2,1,2];plot(angle(H)); axis([−5,140,−4,4]); S_T-sprintf(‘Phase of the impulsereponse function’); h_title=title[S_T];set(h_title,‘Color’,‘k’,‘FontName’,‘Times’,‘FontSize’,10); %d) Mark thevalues that were used by the FPT method hold on; for k=1:r1  x1=p_need(k);   h_o=plot(x1,angle(H(x1)),‘o’);   set(h_o,‘Color’,[1 00]); end; % print -dps filterPT_2.ps % Print separately the original andresults of realization the % linear convolution by the fast paired andrest methods. f_dpt3-figure(‘Position’,[10 398 560 420], ...  ‘Tag’,‘FigD’, ...   ‘Name’,‘ 1-D Signal decomposition v.1, Art M.G.,2002’); h_1=subplot[3,1,1]; plot|f}; axis|[0,260,160,185]];text{230,158,‘t’|; h_title-title(‘Original signal of length 256’|;set(h_title,‘FontName’,‘Times’,‘FontSize’,9); set(h_1,‘FontSize’,9);h_2=subplot(3,1,Z|; plot(fftshift[f_now]}; axis(|0,260,160,185]};S_T=sprintf[(‘The signal filtered by the DPT’); h_title=title[S_T];set(h_title,‘Color’,‘k’,‘FontName’,‘Times’,‘FontSize’,9);set(h_2,‘FontSize’,9); h_3=subplot (3,1,3); plot[fe_real];axis([0,260,160,185]}; S_T=sprintf(‘The signal filtered by the DFT’);h_title=title[S_T]; set(h_title,‘Color’,‘k’,‘FontName’,‘Times’,‘FontSize’,9);set(h_3,‘FontSize’,9); % print -dps filterPT_3.pa

Example 5

Shown here is an example of a MatLab program which implements apair-transform calculation and an inverse transform calculation, inaccordance with embodiments of this disclosure.

% Call; B=fpt_1d(A) % 1-D fast direct paired transform % A --> B oflength N=2**M, M>1. % % Artyom Grigoryan,Yetevan 1996  functionB=fpt_1d(A);  ND=length(A);  ND=log2(ND);  B=A;  LK=D;  NK=ND;  I=1; 1I=1;  while I <= MD   NK=NK/2;   LK−LK+NK;   J=II;   while J <- LK   J1=J+NK;    T−B|J1);    T1=B(J);    B(J1)=T1+T;    B(J)=T1−T;   J=J+1;   end;   II=LK+1;   I−I+1; end; % Call; A=fpti_1d(B) % 1-Dfast inverse paired transform % B --> A of length N=2**M, M>1 % % ArtyumGrigoryan,  Yerevan 1996  function A=fpti_1d[B];   A=B;   ND=length[A];  MD=log2 (ND);   II=ND−1; NK−1;   I−1;   while 1 <= MD     LK=ND+NK;J=II;     while J <- LK       J1=J+NK;       T1=A(J);       T2=A(J1);      T=T1+T2;       T1=T2−T1;       A(J)=T/2.;       A(J1)=T1/2.;      J=J+1;     end;     NK=NK+NK;     II=II−NK;     I=I+1;   end;

REFERENCES

Each of the following references is hereby incorporated in theirentirety.

-   1. A. V. Oppenheim and R. W. Shafer, Digital Signal Processing,    Prentice-Hall, Englewood Cliffs, N.J.: Prentice-Hall, 1989.-   2. J. G. Prookis and D. G. Manolakis, Digital Signal Processing:    Principles, Algorithms, and Applications, 3^(rd) ed., Englewood    Cliffs, N.J.: Prentice-Hall, 1996-   3. O. Ersoy, Fourier-Related Transforms, Fast Algorithms and    Applications, Englewood Cliffs, N.J.: Prentice-Hall, 1997.-   4. C. S. Burrus and T. W. Parks, DFT/FFT and Convolution Algorithms:    Theory and Implementation, New York: Wiley, 1985.-   5. M. T. Heidemann, Multiplicative Complexity, Convolution, and the    DFT, New York: Springer-Verlag, 1988.-   6. A. V. Oppenheim A. V. and R. W. Shafer, Digital signal    processing, Prentice-Hall, Englewood Cliffs, N.J., 1989.-   7. J. G. Prookis and D. G. Manolakis, Digital signal processing:    Principles, Algorithms, and Applications, 3rd ed., Prentice-Hall    Inc., 1996.-   8. O. Ersoy, Fourier-Related Transform, Fast Algorithms and    Applications. Prentice Hall, 1997.-   9. C. S. Burrus and T. W. Parks, DFT/FFT and convolution algorithms:    Theory and Implementation, New York, Wiley 1985.-   10. M. T. Heidemann, Multiplicative complexity, convolution, and the    DFT, New York: Springer-Verlag, 1988.-   11. H. J. Nussbaumer, Fast Fourier transform and convolution    algorithms, 2nd ed., Springer-Verlag, Berlin, 1982.-   12. J. G. Prookis and D. G. Manolakis, Digital signal processing:    Principles, Algorithms, and Applications, 3rd ed., Prentice-Hall    Inc., 1996.-   13. T. G. Stockham, Highspeed convolution and correlation, in 1966    Spring Joint Computer Conf., AFIPS Proc. 28, pp. 229-233.-   14. B. Gold, C. M. Rader, A. V. Oppenheim, T. G. Stockham, Digital    processing of signals, McGraw-Hill, New York 1969, pp. 203-213-   15. R. C. Agarwal, J. W. Cooley, New algorithms for digital    convolution, in 1977 Intern. Conf., Acoust., Speech, Signal    Processing Proc., p. 360.-   16. R. C. Agarwal, C. S. Burrus, Fast one-dimensional digital    convolution by multidimensional techniques, IEEE Trans. ASSP-22, pp.    1-10, 1974.-   17. A. M. Grigoryan, “New algorithms for calculating the discrete    Fourier transforms,” Vichislit. Matem. i Mat. Fiziki, AS USSR, 1986,    vol. 26, no. 5, pp. 84-88.-   18. A. M. Grigoryan, “2-D and 1-D multipaired transforms:    Frequency-time type wavelets,” IEEE Trans. Signal Processing, vol.    49, no. 2, pp. 344-353, February 2001.-   19. A. M. Grigoryan and S. S. Agaian, “Split manageable efficient    algorithm for Fourier and Hadamard transforms,” IEEE Transaction on    Signal Processing, January 2000, vol. 72, no. 1, pp. 172-183.-   20. A. M. Grigoryan, “New algorithms for calculating the discrete    Fourier transforms,” Vichislit. Matem. i Mat. Fiziki, AS USSR, 1986,    vol. 26, no. 5, pp. 84-88.

1. A method for cyclic convolution comprising: obtaining a signal;calculating a paired transform of the signal; grouping components of thepaired transform to form a plurality of splitting-signals, wherein eachof the splitting-signals comprises a plurality of elements; shifting theplurality of elements of each of the plurality of splitting signals tocreate a plurality of shifted signals; multiplying the plurality ofshifted signals by a plurality of moduli of corresponding Fouriertransforms to create a result; and calculating an inverse pairedtransform of the result to filter the signal.
 2. The method of claim 1,wherein the paired transform is calculated in accordance with equation:{{χ′₂ ^(n) _(,2) ^(n) _(t) ; n=0:(r−1), t=0:(2^(r−n−1)−1)}, 1} where nequals log₂ (N), N is a length of the signal, t is an integer variablethat runs from 0 to N−1, r is another integer variable that runs from 0to n−1, and χ′_(p,t) is determined in accordance with equation:χ′_(p,t)(n)=χ_(p,t)(n)−χ_(p,t+N/2)(n) where χ_(p,t) is characterized as:${x_{p,t}(n)} = \{ {{\begin{matrix}{1,} & {{{{if}\mspace{14mu}{np}} = {t\mspace{14mu}{mod}\mspace{14mu} N}},} \\{0,} & {{otherwise},}\end{matrix}\mspace{31mu} n} = {0:{( {N - 1} ).}}} $3. The method of claim 1, wherein the plurality of splitting signals aredetermined for points p from 0 to N−1, where N is a length of thesignal.
 4. The method of claim 3, wherein points p are related to r inaccordance with the equation:p=2^(r) , r=0:(n−1), (and p=0) where n is equal to log₂ N.
 5. The methodof claim 1, wherein each of the plurality of splitting signals areshifted by φ elements, where φ is a value of a normalized phase of oneof the plurality of Fourier transforms.
 6. The method of claim 5,wherein the value of the normalized phase is determined in accordancewith equation:φ=φ/(2π/N) where φ is a phase of the Fourier transform, and N is alength of the signal.
 7. The method of claim 5, wherein shifting furthercomprises shifting the plurality of elements to the right when φ isnegative, and shifting the plurality of elements to the left when φ ispositive.
 8. The method of claim 1, wherein shifting further comprisesreversing the sign of at least one of the plurality of elements.
 9. Themethod of claim 1, wherein calculating an inverse paired transformcomprises calculating an inverse to an unitary paired transforms. 10.The method of claim 9, further comprising multiplying the inverseunitary paired transform by the result.
 11. A computer readable mediacomprising instructions encoded therein for causing a computer system toperform steps comprising: calculating a paired transform of a signal;grouping components of the paired transform to form a plurality ofsplitting-signals; shifting the plurality of splitting signals;multiplying the plurality of splitting signals by a plurality ofcorresponding Fourier transforms; and calculating an inverse pairedtransform of the plurality of splitting signals.
 12. A method for cyclicconvolution comprising: obtaining a signal; steps for calculating apaired transform of the signal; steps for grouping components of thepaired transform to form a plurality of splitting-signals; steps forshifting the plurality of splitting signals; steps for multiplying theplurality of splitting signals by a plurality of moduli of correspondingFourier transforms; and steps for calculating an inverse pairedtransform of the plurality of splitting signals to filter the signal.